function Mat= vdm21(d)
m = (d+1)*(d+2)/2;
[I,J,K] = indices(d);
IM = diag(I)*ones(m,m);
JM = diag(J)*ones(m,m);
KM = diag(K)*ones(m,m);
Mat = (IM/d).^(IM').*(JM/d).^(JM').*(KM/d).^(KM');
IF = gamma(I+1);
JF = gamma(J+1);
KF = gamma(K+1);
A = factorial(d)*ones(m,m)*diag(1./(IF.*JF.*KF));
Mat = A.*Mat;